The present application includes a compact disc containing Appendix A, which is one example of a source code routine written in the C programming language for performing the Gabor spectrogam according to the present invention. The compact disc contains a single text file calle xe2x80x9cGabor Spectrogramxe2x80x9d. The compact disc is in the IBM-PC format and may be viewed with the Windows operating system. The source code listing in the file xe2x80x9cGabor Spectrogramxe2x80x9d is called xe2x80x9cFast Gabor Spectrogram function, Gaborspectrogram.cxe2x80x9d, has a size of 74 kbytes, and a date of creation of Mar. 5, 1999. The source code listing on the compact disc filed in the application herewith is hereby incorporated by reference as though fully and completely set forth herein.
A portion of the disclosure of this patent document, specifically Appendix A, contains material to which the claim of copyright protection is made. The copyright owner has no objection to the facsimile reproduction by any person of the patent document or the patent disclosure, as it appears in the U.S. Patent and Trademark Office file or records, but reserves all other rights whatsoever.
1. Field of the Invention
The present invention relates to digital signal processing and joint time frequency analysis; and more particularly, to systems and methods for determining characteristics of signals having frequency components that vary in time.
2. Description of the Related Art
Due to physical limitations, systems (e.g., an engine, or the human body) or natural phenomena can generally only be studied through the signals (e.g., sound, temperature, heartbeat, and blood pressure) that are generated from or by the system under consideration. Therefore, signal analysis plays a fundamental role in our everyday life. By properly applying analysis techniques, a great deal of information can be obtained about systems without physically invasive procedures.
Most current signal analysis techniques characterize the signal in either the time domain or frequency domain. The time waveform, such as the sound of an engine or an electrocardiogram (ECG), illustrates how the signal""s magnitude varies with time, and the frequency function, such as the power spectrum, indicates how often such changes take place. In most applications, the time and frequency behaviors of signals are closely related to each other. However, conventional techniques are designed to analyze them separately.
Prior art systems for analyzing signals with frequency components that vary in time have used the so-called short time Fourier transform STFT, also known as the windowed Fourier transform. This algorithm is based on a computationally intensive Fourier transform of a large number of short time windows of the input signal. Transforms of each of the short time windows are combined to generate a time varying spectrum of the input signal.
Commonly, the STFT-based spectrogram, described above, has been performed to map time domain (or frequency domain) signals into the joint time-frequency domain. The STFT is simple, but it suffers from the co-called windowing effect, i.e., the STFT is subject to the length of the analysis window function. A short window yields good time resolution, but poor frequency resolution, and vice versa. Particularly, the STFT-based spectrogram is not suitable for the instantaneous frequency estimation which is one of the most important signal""s aspects of interest. The instantaneous frequency estimation from the STFT-based spectrogram is biased, which is subject to the window selection. Finally, the algorithm of the STFT-based spectrogram is not convenient to zoom in to the time-frequency region of interest.
As an alternative, U.S. Pat. No. 5,353,233 to Qian, et al. titled xe2x80x9cMethod and Apparatus for Time Varying Spectrum Analysisxe2x80x9d discloses a technique which may be referred to as a Gabor spectrogram. The Gabor spectrogram is a Gabor expansion based spectrogram. See, Qian, et al., xe2x80x9cDiscrete Gabor Transformxe2x80x9d, IEEE Trans. Signal Processing, Vol.41, No.7, July 1993, pp.2429-2439. Compared to the STFT-based spectrogram, the Gabor spectrogram has better time-frequency resolution and good zoom in capability. The Gabor spectrogram is extremely powerful when the high-resolution time-dependent spectrum is a must. FIG. 1 illustrates time, frequency, and joint time-frequency (Gabor Spectrogram) plots of the sound associated with an aneurysm in the human body. It has been determined that such sound is due to vibration stimulated by the blood flow inside the aneurysm and nearby blood vessels. The mechanism that emits such sounds is complicated and involves the wall, chamber, surrounding blood vessels, and moving blood, under varying pressure. The sound signal is non-stationary, i.e., its frequency contents change with time. Moreover, such a sound record is usually combined with the biological noises generated by the heart, respiratory system, and even eye movement. The bottom plot shows a time waveform that was directly recorded from an intracranial aneurysm during surgery. As expected, the time waveform at the bottom of the figure is rather noisy and thereby virtually provides no useful information for the diagnostic. The spectrum of the signal is illustrated to the left of the figure. The frequency spectrum shows the range of resonant frequencies (500 to 620 Hz), but does not give valuable diagnostic information. When the time domain sound wave is converted into the joint time-frequency domain (the middle plot), it can be clearly seen how the power spectrum evolves over time. Researchers have found that it is much easier to identify the existence of the aneurysm from the joint time-frequency domain than from the time or frequency domain alone.
However, computing the Gabor spectrogram usually takes a long time. Therefore, improved methods are desired for more efficiently or more quickly computing the Gabor spectrogram.
Theoretical Background of the Gabor Spectrogram
The Gabor spectrogram is developed from the decomposition of the Wigner-Ville distribution (WVD). For a given signal s(t), the corresponding Gabor expansion is                               s          ⁡                      (            t            )                          =                              ∑                          m              =                              -                ∞                                      ∞                    ⁢                      xe2x80x83                    ⁢                                    ∑                              n                =                                  -                  ∞                                            ∞                        ⁢                          xe2x80x83                        ⁢                                          C                                  m                  ,                  n                                            ⁢                                                h                                      m                    ,                    n                                                  ⁡                                  (                  t                  )                                                                                        (        1        )            
where Cm,n are known as the Gabor coefficients. The Gabor expansion maps the time domain signal s(t) into a two-dimensional lattice Cm,n. The elementary function hm,n(t) is a time- and frequency-shifted Gaussian function, i.e.,
hm,n(t)=g(txe2x88x92mT)exe2x88x92jnxcexa9txe2x80x83xe2x80x83(2)
where the normalized Gaussian function g(t) has a form                               g          ⁡                      (            t            )                          =                                            (                              πσ                2                            )                                      -              0.25                                ⁢          exp          ⁢                      {                                          -                                  1                                      2                    ⁢                                          σ                      2                                                                                  ⁢                              t                2                                      }                                              (        3        )            
Taking the Wigner-Ville distribution of both sides of Eq.(1) yields                                           WVD            s                    ⁡                      (                          t              ,              ω                        )                          =                              ∑                          m              =                              -                ∞                                      ∞                    ⁢                      xe2x80x83                    ⁢                                    ∑                                                xe2x80x83                                ⁢                                                      m                    xe2x80x2                                    =                                      -                    ∞                                                              ∞                        ⁢                          xe2x80x83                        ⁢                                          ∑                                  n                  =                                      -                    ∞                                                  ∞                            ⁢                              xe2x80x83                            ⁢                                                ∑                                                            n                      xe2x80x2                                        =                                          -                      ∞                                                        ∞                                ⁢                                  xe2x80x83                                ⁢                                                      C                                          m                      ,                      n                                                        ⁢                                      xe2x80x83                                    ⁢                                      C                                                                  m                        xe2x80x2                                            ,                                              n                        xe2x80x2                                                              *                                    ⁢                                                            WVD                                              h                        ,                                                  h                          xe2x80x2                                                                                      ⁡                                          (                                              t                        ,                        ω                                            )                                                                                                                              (        4        )            
where WVDh,hxe2x80x2(t,xcfx89), the WVD of the time- and frequency-shifted Gaussian function, which can be referred to as an xe2x80x9cenergy atomxe2x80x9d, is concentrated, oscillated, and symmetrical. The energy atom has a closed form, and thus can be pre-calculated and stored in a table. If the Wigner-Ville distribution WVDs(t,xcfx89) describes the signal""s energy distribution in the joint time-frequency domain, then Eq.(4) shows that the signal energy is a superposition of an infinite number of energy atoms. Note that the energy (average) contained in each individual WVDh,hxe2x80x2(t,xcfx89) is inversely proportional to the distance between corresponding elementary functions hm,n(t) and hmxe2x80x2,nxe2x80x2(t). Based on their contribution to the entire signal energy, Eq.(4) can be re-grouped to obtain the Gabor spectrogram as                                           GS            D                    ⁡                      (                          t              ,              ω                        )                          =                              ∑                                                            "LeftBracketingBar"                                      m                    -                                          m                      xe2x80x2                                                        "RightBracketingBar"                                +                                  "LeftBracketingBar"                                      n                    -                                          n                      xe2x80x2                                                        "RightBracketingBar"                                            ≤              D                                ⁢                      xe2x80x83                    ⁢                                    C                              m                ,                n                                      ⁢                          C                                                m                  xe2x80x2                                ,                                  n                  xe2x80x2                                            *                        ⁢                                          WVD                                  h                  ,                                      h                    xe2x80x2                                                              ⁡                              (                                  t                  ,                  ω                                )                                                                        (        5        )            
The parameter D denotes the order of the Gabor spectrogram, which characterizes the Manhattan distance between two elementary functions, hm,n(t) and hmxe2x80x2,nxe2x80x2(t). The farther the elementary functions, hm,n(t) and hmxe2x80x2,nxe2x80x2(t), are from each other, the stronger the energy atom WVDh,hxe2x80x2(t,xcfx89) oscillates.
For the discrete version, the time and frequency indices of GSD(t,xcfx89) in Eq.(5) can be quantized, i.e.,                                           GS            D                    ⁡                      [                          i              ,              k                        ]                          =                              ∑                                                            "LeftBracketingBar"                                      m                    -                                          m                      xe2x80x2                                                        "RightBracketingBar"                                +                                  "LeftBracketingBar"                                      n                    -                                          n                      xe2x80x2                                                        "RightBracketingBar"                                            ≤              D                                ⁢                      xe2x80x83                    ⁢                                    C                              m                ,                n                                      ⁢                          C                                                m                  xe2x80x2                                ,                                  n                  xe2x80x2                                            *                        ⁢                                          WVD                                  h                  ,                                      h                    xe2x80x2                                                              ⁡                              [                                                      i                    ⁢                                          xe2x80x83                                        ⁢                                          Δ                      t                                                        ,                                      k                    ⁢                                          xe2x80x83                                        ⁢                                          Δ                      f                                                                      ]                                                                        (        6        )            
where xcex94t and xcex94f denote the time and frequency sampling intervals, respectively. By varying sampling intervals, zoom in and zoom out can be easily performed. This is one of the major advantages of the Gabor spectrogram over other joint time-frequency analysis techniques, such as the STFT and the Wigner-Ville distribution. For the STFT or the WVD, special filters must be applied to zoom in or out. Besides the difficulty of the filter design, the main disadvantage of using filters is that the resulting presentation may not be consistent. For example, the quantity of STFT or WVD at the same time and frequency point could vary at different scales due to the different filters involved.
The current algorithm first computes the Gabor coefficients Cm,n and then builds the Gabor spectrogram simply by looking at values in the table. The speed and accuracy of this approach relies heavily on the size of the table (pre-calculated energy atoms). The bigger the table, the greater the accuracy. On the other hand, however, the bigger the table, the more elements of the matrix GSD[i,k] that need to be updated. In other words, a larger table (more accuracy) requires longer processing time. For the sake of computation speed, the precision of the table of energy atoms is usually limited to be no less than 10xe2x88x923, which has been found to be insufficient for many applications.
Therefore, an improved system and method is desired for more efficiently, more quickly, and/or more accurately computing the Gabor spectrogram.
The present invention provides a signal analyzer, method and memory medium for generating a time varying spectrum for input signals characterized by time-varying frequencies. The invention can also be characterized as a method for analyzing a signal based on the time varying spectrum. The signal analyzer includes a source of a sequence of digital signals representative of an input signal, a processor coupled to the source, and a memory medium coupled to the processor. The memory medium stores a software program which is executable by the processor to compute the time-varying spectrum of the input signal.
When the processor executes the software program, the processor is operable to first compute a Gabor transform (that is, a sampled short-time Founrer transform) of the digital signals to produce Gabor coefficients. The processor then computes an auto-correlation of the Gabor coefficients to produce auto-correlation results. The auto-correlation results are then applied to a 2-dimensional interpolation filter to produce the time-varying spectrum, wherein the time-varying spectrum is a Gabor spectrogram. The signal analyzer may repeat the above steps n+1 times and sum the results for an n order time-varying spectrum. The process may then operate to process and/or display the time-varying spectrum.
In computing the Gabor coefficients of the digital signals, the processor is operable to compute orthogonal-like discrete Gabor transform Cm,n in response to the digital signals. The processor computes the Gabor transform of the digital signals by computing:                                                                         C                                  m                  ,                  n                                            =                              xe2x80x83                            ⁢                                                ∑                  k                                ⁢                                  xe2x80x83                                ⁢                                                      s                    ⁡                                          [                                              k                        +                        mT                                            ]                                                        ⁢                                      γ                    ⁡                                          [                      k                      ]                                                        ⁢                                      W                                          L                      /                      Ω                                                              -                                              n                        ⁡                                                  [                                                      k                            +                            mT                                                    ]                                                                                                                                                                                            =                              xe2x80x83                            ⁢                                                W                  N                                      -                    nmT                                                  ⁢                                                      ∑                                          l                      =                      0                                                                                      L                        /                        N                                            -                      1                                                        ⁢                                      xe2x80x83                                    ⁢                                                            ∑                                              k                        =                        0                                                                    N                        -                        1                                                              ⁢                                          xe2x80x83                                        ⁢                                                                  s                        ⁡                                                  [                                                      lN                            +                            k                            +                            mT                                                    ]                                                                    ⁢                                              γ                        ⁡                                                  [                                                      lN                            +                            k                                                    ]                                                                    ⁢                                              W                        N                                                  -                          nk                                                                                                                                                                                            =                              xe2x80x83                            ⁢                                                W                                      N                    /                    T                                                        -                    nm                                                  ⁢                                                      ∑                                          l                      =                      0                                                                                      L                        /                        N                                            -                      1                                                        ⁢                                      xe2x80x83                                    ⁢                                                            ∑                                              k                        =                        0                                                                    L                        -                        1                                                              ⁢                                          xe2x80x83                                        ⁢                                                                                            ℜ                                                      l                            ,                            m                                                                          ⁡                                                  [                          k                          ]                                                                    ⁢                                              W                        N                                                  -                          nk                                                                                                                                                                            (        18        )            
wherein N denotes the number of frequency bins,
wherein L is the length of the filter, which is greater than or equal to N.
In applying the auto-correlation results to the interpolation filter to produce the time-varying spectrum, the processor is operable to compute:                               y          ⁡                      [            n            ]                          =                              ∑                          m              =              0                        M                    ⁢                      xe2x80x83                    ⁢                                    x              ⁡                              [                m                ]                                      ⁢                          h              ⁡                              [                                                                            (                                              n                        -                                                  n                          0                                                                    )                                        ⁢                    Δ                                    -                  mT                                ]                                                                        (        12        )            
where the constants n0, xcex94, and T denote time delay, decimation, and interpolation factors. The memory medium may store a table of pre-computed filter response values, and these pre-computed filter response values may be used in applying the auto-correlation results to the interpolation filter to produce the time-varying spectrum. The decimation and interpolation filter described in (12) can be efficiently computed in the frequency domain, i.e.,                               Y          ⁡                      (            θ            )                          =                              ⅇ                                          -                                                      ⅈ                    ⁢                    n                                    0                                            ⁢              θ                                ⁢                      X            ⁡                          (                                                T                  Δ                                ⁢                θ                            )                                ⁢                      H            ⁡                          (                              θ                Δ                            )                                                          (        13        )            
where Y(xcex8) denotes the Fourier transform of y[n].
The previous algorithm of the Gabor spectrogram usually employs a look-up table. The method presented in this patent converts that process into the standard two dimensional multi-rate filter operation. The resulting method is faster and amenable for hardware implementation.